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We investigate the domain structure of pair contact process with diffusion (PCPD). PCPD is a 
stochastic reaction-diffusion model which evolves by the competition of two binary reactions, 2A — » 
3A and 2A — > 0. In addition, each particle diffuses isotropically, which leads to the bidirectional 
coupling between solitary particles and pairs. The coupling from pairs to solitary particles is linear, 
while the opposite coupling is quadratic. The spreading domain formed from localized activities 
in vacuum consists of two regions, the coupled region of size R p where pairs and solitary particles 
coexist and the uncoupled region of size Ru where only solitary particles exist respectively. As the 
size of the whole domain R is given as R — R p + Ru, R P and Ru are the basic length scales of 
PCPD. At criticality, R p and Ru scale as R p ~ t 1/Zp and Ru ~ t 1/Zu with Z v > Z p . We estimate 
Z p = 1.61(f) and Zu = f. 768(8). Hence, the correction to the scaling of R, Q — Ru/R P extremely 
slowly decays, which makes it practically impossible to identify the asymptotic scaling behavior of 
R. In addition to the generic feature of the bidirectional coupling, the double domain structure is 
another reason for the extremely slow approach to the asymptotic scaling regime of PCPD. 

PACS numbers: 64.60.-i, 05.40.-a, 82.20.-w, 05.70.Ln 



I. INTRODUCTION 

The concept of universality makes possible to un- 
derstand and classify the various and complicate crit- 
ical behavior of a number of equilibrium models 
For nonequilibrium critical phenomena, extensive studies 
during past decades revealed that various nonequilibrium 
systems also exhibit universal behavior characterized by 
several features |2| . Hence, it has been an important issue 
to identify nonequilibrim universality classes by finding 
common physical features. 

Among nonequilibrium critical phenomena, absorbing 
phase transitions (APT's) from fluctuating active states 
into absorbing states in which the system is trapped for- 
ever have been a field of growing interest during last 
decades 0, OJ. Recent theoretical and numerical stud- 
ies show that APT's exhibit universality and it can be 
classified according to conservation laws, dimensionality 
of systems and symmetries of absorbing states 0, H, [j]. 
However only a few universality classes have been iden- 
tified so far. Directed percolation (DP) Q|, [|, JE [jf and 
parity conserving (PC) 0, H, H EE EH E2> EH class are 
well studied classes among others. DP class includes sys- 
tems with no special attributes except the time reversal 
symmetry, so that most systems studied so far belong to 
this class. 

As a research direction to search for further unknown 
universality classes, coupled s yst ems have been studied 

recently @ Q EE til ES El, M M M Hi HI HI HE 

. A coupled system is a multi-species system in which 
each species is coupled to the others in certain ways such 
as bidirectional and unidirectional coupling in linear or 
quadratic ways. However the coupled systems do not al- 
ways exhibit new critical behavior. For bidirectionally 
coupled systems, the critical behavior depends on the 
manner of the coupling. For instance, quadratically cou- 
pled DP systems still belong to DP class despite their 



complex behavior [14| . However, linearly coupled sys- 
tems belonging to DP or PC class exhibit mean-field or 
non-trivial critical behavior [9|, Il5l . Il6l . |17L |18| | . Linearly 
and unidirectionally coupled systems exhibit new critical 
behavior at multicritical point where all sub-systems are 

critical 0, HE HH M, HE HI HE Hi . 

Among single species systems, pair contact process 
with diffusion (PCPD) can be regarded as a two species 
system. PCPD has been extensively studied during last 
years due to its nontrivial critical behavior (see [13] for re- 
view). However, in spite of extensive theoretical and nu- 
merical studies, the critical behavior is not clearly uncov- 
ered yet. PCPD is a stochastic reaction-diffusion model, 
which evolves by the competition of two processes, fis- 
sion (2 A — > 3 A) and annihilation (2 A — > 0). In addi- 
tion, each particle performs isotropic diffusion. With- 
out diffusion, the model is so-called pair contact process 
(PCP) belonging to DP class (28[. Since the reactions 
involve pairs, diffusing solitary particles are not engaged 
in the binary reactions. However when two solitary par- 
ticles form a pair, the reactions take place. On the other 
hand, solitary particles are created from pairs by diffu- 
sion. Hence, PCPD can be regarded as a bidirection- 
ally coupled two species system in which the order of the 
coupling is linear in the direction from pairs to solitary 
particles and quadratic in the opposite direction. This 
observation leads to the cyclically coupled DP and pair 
annihilation which exhibits the similar type of critical 
behavior to that of PCPD 

In this paper, we investigate the domain structure of 
PCPD. When a spreading domain is formed from local- 
ized initial activities, the quadratic coupling from solitary 
particles to pairs allows the pair-free region in which only 
solitary particles are present. We call the pair-free region 
so-called uncoupled region. On the other hand, the lin- 
ear coupling from pairs to solitary particles results in the 
so-called coupled region in which pairs and solitary par- 
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FIG. 1: (Color online) The double domain structure of PCPD. 
The region with size R p is the coupled region in which pairs 
(blue) and solitary particles (red) coexist. The region with 
size Ru is the uncoupled region in which only solitary particles 
exist. 



tides coexist. The pair- free region encloses the coupled 
region as shown in Fig. [TJ Hence, the spreading domain 
is divided into two regions, coupled and uncoupled re- 
gion. This kind of double domain structure was found 
in unidirectionally coupled two level hierarchies [IE H|| . 
As shown in previous studies on unidirectionally coupled 
systems [H, EE EH, HE Efi| , the measurements of criti- 
cal exponents are very difficult due to long-time drift of 
the exponents. The one reason of the drift is the generic 
feature of the unidirectional couplin g [22I |23| , the other 
is the double domain structure [25L l26j| . In measuring 
critical exponents, one can overcome the latter effect by 
measuring quantities in each region separately. 

To explain the effect of the domain structure on the 
critical behavior, let us consider the unidirectionally cou- 
pled DP-PC coupling HU. In the DP-PC coupling, PC 
process (B) is unidirectionally coupled with DP process 
(A) via A — > A + 2B (linear). Since DP process spreads 
faster than PC process, the coupling changes the spread- 
ing behavior of PC process. There are no interactions 
between A and B species. Starting with a single A par- 
ticle on A-level, a P-domain is form and spreads faster 
than a A-domain does due to the coupling. As a result, 
the uncoupled region where no A particles exist is formed 
[25l ]. The coupled region is just the region of P-domain 
overlapped with the >l-domain where the coupling ex- 
ists. Hence, the size of the P-domain (Rb) is the sum 
of the size of coupled (Rc) and uncoupled region (Ru), 
Rb = Rc + Ru- At the multicritical point where both 
species are critical, the size of each domain increases in 
power-law such as Ra ~ t ZA , Rb ~ t ZB , Rc ~ t zc and 
Ru ~ t Zu . Since zc — za by definition, the scaling of 
Rb is given as Rb ~ at ZA + bt Zu , a and b are constants. 
A recent study showed za > Zu for the fast spreading A- 
domain (source) [2^|. Hence, Rb scales as t ZA asymptoti- 
cally for faster spreading source A, and Ru plays the cor- 
rection to the scaling of Rb- Since the difference za — zu 
is sufficiently small in the DP-PC coupling, the existence 
of the uncoupled region makes precise estimates of zb 
difficult within moderate simulation time. For the DP- 
PC coupling, the estimate of 2zb is 1.24(1) which is less 
than the expectation 2zb = 2zdp — 1.265 (2(|. On the 
other hand, the estimate of 2zu is 2zu — 1.20(2) which 
is compatible with the leading scaling. 

As in unidirectionally coupled systems, it is expected 



that the double domain structure of PCPD also makes 
it difficult to identify the critical spreading behavior pre- 
cisely The aim of this paper is to investigate the effect of 
the double domain structure on the critical spreading of 
PCPD and to estimate critical exponents more precisely. 
We consider the PCPD of Ref. [30]. In this model, a 
randomly selected particle attempts to hop to one of the 
nearest neighbor sites with an equal probability. If the 
target site is empty, the attempt is accepted. However 
if the target site is occupied, (i) two particles annihilate 
with probability p or (ii) the hopping attempt is rejected 
and the pair (the chosen particle and one at the target 
site) tries to create a particle at the randomly chosen 
nearest neighbor site of the pair. When the target site is 
occupied, the branching attempt is rejected. The critical 
point of this model is p c = 0.133 519(3) in one dimension 

PCPD has three sectors in configuration spaces accord- 
ing to the existence of pairs(P) and solitary particles(S). 
The one is the configurations in which both pairs and soli- 
tary particles are present (PS-ensemble). We call con- 
figurations with at least one pair (two solitary particles) 
P-ensemble (S-ensemble). In P-ensemble (S-ensemble), 
solitary particles (pairs) may be present or not. The P- 
ensemble is the reactive subspace of Ref. [3l| . A conven- 
tional ensemble includes configurations with at least two 
particles which can be either two solitary particles or one 
pair. We call the conventional ensemble All-ensemb\e. 
Since solitary particles are effectively linearly coupled to 
pairs, the existence of pairs implies the existence of soli- 
tary particles. Hence, PS-ensemble should coincide with 
P-ensemble asymptotically. However, solitary particles 
transform into pairs by collisions, the coupling in this di- 
rection is quadratic. Hence, the existence of solitary par- 
ticles does not always guarantee the presence of pairs due 
to the long life time of solitary particles. So S'-ensemble 
coincide with All-ensemb\e. As a result, there are two 
distinct ensembles in PCPD. 

We define the size of each domain as follows (See 
Fig. [T|). The size of the whole domain at time t (R(t)) 
is defined as the distance between the leftmost and the 
rightmost particle. When both pairs and solitary parti- 
cles exist simultaneously, we can define the size of the 
coupled region (Rc(t)) and the size of the uncoupled re- 
gion (Ru(t j). Since solitary particles are linearly coupled 
to pairs, we define the size of the coupled region (Rc(t)) 
as the spreading distance of pairs (R p ) defined as the 
distance between the leftmost and the rightmost pair. 
Then, Ru(t) is given as Ru = R — Rp. Hence, we have 
three different lengths, R, R p and Ru in PCPD. To take 
into account the three length scales at the same time, one 
should use PS'-ensemble in which only two lengths, R p 
and Ru are the fundamental length scales of PCPD due 
to R = Rp + Ru- 

At p c — 0.133 519, we perform defect Monte Carlo 
simulations with a pair on one dimensional empty lat- 
tice. We run simulations up to t = 10 time steps using 
3.6 x 10 6 independent runs. We measure the squared 
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FIG. 2: (Color online) The scaling plots of various sizes, (a) 
The scaling plot of R 2 (black) and R v (blue) with Z = 1.663 
and Zjj = 1.768. (b) the scaling plot of R p with Z v = 1.61. 




FIG. 3: The scaling plot of the ratio Q = Ru /Rp- The main 
plot shows Qt with (f> = 0.054. The inset shows the double 
logarithmic plot of Q(t). 



sizes, R 2 , R 2 and R\j for surviving PS-ensemble. The 
PCPD has two absorbing states, vacuum and states with 
one diffusing solitary particle. Hence, we stop the sim- 
ulations when the total number of particles N is less 
then two. At the criticality, the squared sizes scale 
as R 2 ~ i 2 / z , R 2 ~ t 2 l z * and R% ~ t 2 / z " . Fig. rj 
shows the scaling plots of the squared sizes, R 2 /t 2 / Zn >. 
We obtain the best scaling plot with Z = 1.663(5), 
Z p = 1.61(1) and Z v = 1.768(8) respectively. The er- 
rors of our estimates should be larger due to the error 
of p c . Within the numerical errors at the criticality, our 
estimate of Z agree with the previous studies [27| , espe- 
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FIG. 4: The scaling plot of N p . The main plot shows N v /f 
with rj = 0.275. The inset shows the semilogarithmic plot of 
the ratio N a /N p . 



daily Z — 1.70(5) [T(|. Since Zu > Z p , the total spread- 
ing distance R(= R p + Ru) should scale as R ~ t l / Zp , 
and Ru plays the role of the correction to the scaling as 
in unidirectionally coupled systems. Hence, we conclude 
Z = Z p — 1.61(1) which is the smallest value among the 
estimates of previous studies |27| . For reference, we also 
measure R 2 using ^4^-ensemble. The difference of All- 
ensemble from PS-ensemble or P-ensemble is that All- 
ensemble includes configurations without pairs. Since 
solitary particles spread diffusively, it is expected that R 
averaged over ^MZ-ensemble scales differently from that of 
PS-ensemble. From the scaling plot of R 2 /t 2 / z , we esti- 
mate Z — 1.676(3) for ^4^-ensemble which is larger than 
that of PS-ensemble (not shown) . The slow spreading of 
the whole domain in yW-ensemble results from the diffu- 
sive motions of solitary particles in configuration without 
pairs. 

Since R scales as R = R p (l + Ru/R P ), the correction 
to the R is Q — Ru/R P which decays as Q ~ t^^ with 
4> = {Zu — Z p ) I Z p Zu ■ To see how the correction decays 
slowly in time, we plot Q in Fig. [U When Q <g; 1, the 
correction is negligible. However, as shown in the inset, 
Q is still comparable to one even at t = 10 7 . We obtain 
the best scaling plot of Qt with <j) = 0.054(4). As Q 
decays with very small (j>, it is practically impossible to 
reach the asymptotic scaling regime of R ~ R p . As a 
result, one should take the double domain structure into 
account for the more precise measurement of the dynamic 
exponent Z. 

In addition to the sizes, we also measure the number 
of pairs (N p ) and of solitary particles (N s ) averaged over 
all samples, and the survival probability of PS-ensemble 
(P ps ) and All-ensemble (P «). As solitary particles are 
linearly coupled to pairs, N s is proportional to N p . At 
criticality, N p scales as N p ~ t v . Fig. [4] shows the scal- 
ing plot N p /t v . We obtain the best scaling plot with 
rj = 0.275(5). The inset show the ratio of N s /N p , which 
converges to one as expected. P ps {t) and P a ii{t) de- 



4 



&h 0.7 

0.6 

- Oh 

^ 0.5 



"3 




- 


a. 


















0.1 


1 1 1 1 f 




10 2 1 3 1 4 1 5 10* 10 

t 




1 1 


1 1 1 



io 2 io 3 io 4 io 5 io 6 io 7 
t 



FIG. 5: (Color online) The scaling plot of P ps and P a ii- The 
main plot shows the scaling plot of P ps t Spa (red) and P a iit Sa 
(black) with 5' ps = S' a — 0.13. We add a constant —0.4 to 
the scaling plot of P a u for the better presentation. The inset 
shows the double logarithmic plots. 

cay in power-law as P ps ~ t~ s p s and P a ii ~ t~ Sa at 
the criticality. Fig. \E\ shows the scaling plots of P x t s * 
and the double logarithmic plots of P ps and P a ii- With 
dps = fi'a = 0.130(3), we obtain the best scaling plots. 
The system escapes from the non-reactive sector in which 
only diffusing solitary particles are present via pair anni- 
hilations of solitary particles. As an naive argument for 
the equality of S' ps = 5' a , the escaping probability from 

non-reactive sector within time r may scale as 1 — r -1 ' 2 
which is the death probability of two diffusing particles 
undergoing the reaction S + S — > within r. Hence, the 
contribution of non-reactive sector to P a u is negligible 
due to the fast escaping probability, which leads to the 
same scaling behavior of P a ii as P p . 

In summary, we investigate the domain structure of 
PCPD. We numerically confirm the existence of the dou- 
ble domain structure in PCPD. This double domain 
structure comes from the linear and quadratic bidirec- 
tional couplings. The structure intrinsically makes the 
serious correction to the scaling of the critical spreading 
of a domain. 

Starting with a pair, a domain grows and spreads in 
vacuum via fission (A + A — > 3 A) and spontaneous anni- 
hilation of pairs (2 A — ► 0) in PCPD. In addition to the bi- 
nary reactions, each particle diffuses isotropically, which 
leads to the bidirectional coupling between solitary par- 
ticles and pairs. The coupling from pairs to solitary par- 
ticles is linear, while the opposite coupling is quadratic. 
The difference of the coupling ways results in the double 
domain structure of the whole domain, the coupled and 
the uncoupled region respectively. As a result, the size 
of the whole domain (R) is given as the sum of the size 
of the coupled region (R p ) and of the uncoupled region 



(Rjj)- Hence, R p and Rjj are the basic length scales char- 
acterizing the scaling behavior of the spreading domain 
in PCPD. We numerically find that R p and Rjj scale as 
R p ~ i 1 /^ and R v ~ t 1 ^" with Z v > Z c at criticality. 
Hence, R should asymptotically scale as R ~ t x l z with 
Z = Z p and Rjj plays the role of the correction to the 
scaling. However, the direct measurement of R leads to 
the underestimate of the asymptotic value of Z because 
the correction Q = Ru/R p decays with very small ex- 
ponent. Since it is practically impossible to reach the 
asymptotic scaling region of R ~ t 1 / 2 ', it is important 
to take the domain structure into account in simulations 
for more precise estimate of the dynamic exponent Z of 
PCPD. 

We classify particle configurations into four ensembles, 
which are finally reduced to two distinct ensembles, P- 
ensemble and A^-ensemble respectively. All-ensemble 
includes configurations without pairs, while P-ensemble 
does not. The survival probabilities of two ensembles de- 
cay with the same exponent. However, the whole domain 
appears to spread more slowly in All-ensemble than in 
P-ensemble due to the diffusive motions of solitary par- 
ticles in configurations without pairs. Hence, in addition 
to the domain structure, the diffusive motions of solitary 
particles in All-ensemble raise another correction to the 
scaling of the total spreading distance R which does not 
appear in P-ensemble. 

As the linear-quadratic bidirectional coupling is the 
common feature of various PCPD studied so far, the dou- 
ble domain structure should appear in other PCPD vari- 
ants. Among PCPD variants, we investigate the domain 
structure of the bosonic PCPD with soft-constraint of 
Ref. pr| which is known to exhibit the clear power-law 
decays at the criticality. For this model, we also confirm 
the existence of the double domain structure and the crit- 
ical spreading behavior similar to that of PCPD studied 
in this paper. Hence, the double domain structure is a 
common feature of PCPD variants. 

As shown in recent studies, as PCPD extremely slowly 
approaches to its asymptotic scaling region, it is very dif- 
ficult to identify the critical behavior. The bidirectional 
coupling should be the main reason in itself. On the other 
hand, the double domain structure naturally appears due 
to the linear-quadratic bidirectional couplings, which en- 
hance the long-time drift of dynamic exponent Z. One 
can overcomes the latter effect by considering the scaling 
behavior of sub-domains separately as in unidirectionally 
coupled systems. 
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